/*
 * thermal.c
 *
 *  Created on: Mar 12, 2012
 *      Author: qian
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "thermal.h"
#include "stack.h"
#include "auxiliary.h"
#include "/usr/include/suitesparse/klu.h"

/*allocate memory for thermal data*/
thermal_data_t* init_thermal_data(stk_descr_t *stk)
{
	int num_source_cell = stk->rows*stk->cols*stk->n_layers;
	int total_units=obtain_total_units(stk);
	int total_ch_num;
	thermal_data_t *thermo=(thermal_data_t*)malloc(sizeof(thermal_data_t));

	/*allocate memory for common elements*/
	thermo->gx_BEOL=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gy_BEOL=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gz_BEOL=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gx_source=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gy_source=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gz_source=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gx_bulk=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gy_bulk=(double*)malloc(stk->n_layers*sizeof(double));
	thermo->gz_bulk=(double*)malloc(stk->n_layers*sizeof(double));

	thermo->power_grid=(double*)malloc(num_source_cell*sizeof(double));
	thermo->steady_temp_grid=(double*)malloc(num_source_cell*sizeof(double));
	thermo->steady_temp_block=(double*)malloc(total_units*sizeof(double));

	/*allocate memory depending on cooling method*/
	if(stk->cooling_method==1)
		thermo->pk_thermal=(package_thermal_t *)malloc(sizeof(package_thermal_t));
	else if(stk->cooling_method==2)
	{
		total_ch_num=stk->microchannel->n_channels_per_layer*(stk->n_layers-1);
		thermo->ch_thermal=(channel_thermal_t *)malloc(sizeof(channel_thermal_t));
		thermo->ch_thermal->h_ch=(double*)malloc(total_ch_num*stk->rows*sizeof(double));
		thermo->ch_thermal->flowrate=(double*)malloc(total_ch_num*sizeof(double));
		thermo->ch_thermal->Tch=(double*)malloc(total_ch_num*stk->rows*sizeof(double));
		thermo->ch_thermal->ch_gz=(double*)malloc(total_ch_num*stk->rows*sizeof(double));
		thermo->ch_thermal->ch_gx=(double*)malloc(total_ch_num*stk->rows*sizeof(double));
	}
	else
		fatal("Undefined cooling method,please check configuration file!");

	return thermo;
}


/*populate thermal data, fill in power information and produce conductances*/
void fill_thermal_data(stk_descr_t *stk, thermal_data_t *thermo)
{
	/*generate power map*/
	double *power_block;
	power_block=obtain_unit_powers(stk);
	xlate_vector_b2g(stk, thermo, power_block, V_POWER);

	populate_silicon_conductances(stk, thermo);

	if(stk->cooling_method==1)
		populate_package_conductances(stk,thermo);
	else if(stk->cooling_method==2)
		populate_channel_conductances(stk, thermo);
	else
		fatal("undefined cooling method when populate conductances, check config file.");

	thermo->map_mode=GRID_AVG;
}

/*store unit powers into a 1-D array*/
double *obtain_unit_powers(stk_descr_t *stk)
{
	int n_total_units=obtain_total_units(stk);
	int i,j,base;

	double *power;
	power=(double *)malloc(n_total_units*sizeof(double));

	for(i=0,base=0;i<stk->n_layers;i++)
	{
		for(j=0;j<stk->layers[i].flp->n_units;j++)
			power[j+base]=stk->layers[i].flp->units[j].power;
		base+=stk->layers[i].flp->n_units;
	}
	return power;

}

/*fill in conductances regarding die cells*/
void populate_silicon_conductances(stk_descr_t *stk, thermal_data_t *thermo)
{
	int i;
	double BEOL_t, source_t, bulk_t;
	double cell_width=stk->width/stk->cols;
	double cell_length=stk->length/stk->rows;

	for(i=0;i<stk->n_layers;i++)
	{
		BEOL_t=stk->layers[i].BEOL_height;
		source_t=stk->layers[i].source_height;
		bulk_t=stk->layers[i].bulk_height;
		/*gx's and gy's are cell to cell conductances*/
		thermo->gx_BEOL[i]=K_BEOL*cell_length*BEOL_t/cell_width;
		thermo->gy_BEOL[i]=K_BEOL*cell_width*BEOL_t/cell_length;

		thermo->gx_source[i]=K_SILICON*cell_length*source_t/cell_width;
		thermo->gy_source[i]=K_SILICON*cell_width*source_t/cell_length;

		thermo->gx_bulk[i]=K_SILICON*cell_length*bulk_t/cell_width;
		thermo->gy_bulk[i]=K_SILICON*cell_width*bulk_t/cell_length;

		/*gz's are cells' boundary to boundary conductances to cater different cooling architecture*/
		thermo->gz_BEOL[i]=K_BEOL*cell_width*cell_length/BEOL_t;
		thermo->gz_source[i]=K_SILICON*cell_width*cell_length/source_t;
		thermo->gz_bulk[i]=K_SILICON*cell_width*cell_length/bulk_t;
	}

}

/*fill in package conductances*/
void populate_package_conductances(stk_descr_t *stk, thermal_data_t *thermo)
{
	double width = stk->width;
	double height = stk->length;
	double cw=width/stk->cols;
	double ch=height/stk->rows;

	double s_sp = stk->heatsink->width_spreader;
	double t_sp = stk->heatsink->t_spreader;
	double s_sink = stk->heatsink->width_heatsink;
	double t_sink = stk->heatsink->t_heatsink;
	double r_conv = stk->heatsink->r_conv;
	double t_tim = stk->heatsink->t_TIM;

	double k_sink = stk->heatsink->k_heatsink;
	double k_sp = stk->heatsink->k_spreader;
	double k_tim = stk->heatsink->k_TIM;

	/*core part thermal resistances of spreader and sink*/
	thermo->pk_thermal->sp_rx =getr(k_sp, cw, t_sp*ch);
	thermo->pk_thermal->sp_ry =getr(k_sp, ch, t_sp*cw);
	thermo->pk_thermal->sp_rz =getr(k_sp, t_sp/2, cw*ch); //centre to edge

	thermo->pk_thermal->sink_rx = getr(k_sink, cw, t_sink*ch);
	thermo->pk_thermal->sink_ry = getr(k_sink, ch, t_sink*cw);
	thermo->pk_thermal->sink_rz = getr(k_sink, t_sink/2, cw*ch); //centre to edge
	thermo->pk_thermal->sink_ramb=(r_conv*s_sink*s_sink)/(cw*ch);

	/* lateral R's of spreader and sink, all are from centre of the node to its corresponding edge*/
	/* lateral resistances of west and east extra nodes of spreader (to core)*/
	thermo->pk_thermal->rx_sp_per = getr(k_sp, (s_sp-width)/4, (s_sp+3*height)/4 * t_sp);
	/* lateral resistances of north and south extra nodes of spreader (to core)*/
	thermo->pk_thermal->ry_sp_per = getr(k_sp, (s_sp-height)/4, (s_sp+3*width)/4 * t_sp);
	/* lateral resistances of west and east middle extra nodes of heatsink (to core)*/
	thermo->pk_thermal->rx_hs_ci = getr(k_sink, (s_sp-width)/4, (s_sp+3*height)/4 * t_sink);
	/* lateral resistances of north and south middle extra nodes of heatsink (to core)*/
	thermo->pk_thermal->ry_hs_ci = getr(k_sink, (s_sp-height)/4, (s_sp+3*width)/4 * t_sink);
	/* lateral resistances of west and east middle extra nodes of heatsink (to outer extra nodes of heatsink)*/
	thermo->pk_thermal->rx_hs_co = getr(k_sink, (s_sp-width)/4, (3*s_sp+height)/4 * t_sink);
	/* lateral resistances of north and south middle extra nodes of heatsink (to outer extra nodes of heatsink)*/
	thermo->pk_thermal->ry_hs_co = getr(k_sink, (s_sp-height)/4, (3*s_sp+width)/4 * t_sink);
	/* lateral resistances of outer extra nodes of heatsink (to middle extra nodes of heatsink,
	 * four of them are identical since both spreader and heatsink are squarish)*/
	thermo->pk_thermal->rxy_hs_outer = getr(k_sink, (s_sink-s_sp)/4, (s_sink+3*s_sp)/4 * t_sink);

	/* vertical R's of spreader and sink, from centre to edge */
	thermo->pk_thermal->rz_sp_we = getr(k_sp, t_sp/2, (s_sp+height) * (s_sp-width) / 4);
	thermo->pk_thermal->rz_sp_ns = getr(k_sp, t_sp/2, (s_sp+width) * (s_sp-height) / 4);
	thermo->pk_thermal->rz_hs_c_we = getr(k_sink, t_sink/2, (s_sp+height) * (s_sp-width) / 4);
	thermo->pk_thermal->rz_hs_c_ns = getr(k_sink, t_sink/2, (s_sp+width) * (s_sp-height) / 4);
	thermo->pk_thermal->rz_hs_outer = getr(k_sink, t_sink/2, (s_sink*s_sink - s_sp*s_sp) / 4);

	/* vertical R's from heatsink to ambient (divide r_convec proportional to area) */
	thermo->pk_thermal->ramb_c_we = r_conv * (s_sink * s_sink)/((s_sp+height) * (s_sp-width) / 4);
	thermo->pk_thermal->ramb_c_ns = r_conv * (s_sink * s_sink)/((s_sp+width) * (s_sp-height) / 4);
	thermo->pk_thermal->ramb_outer = r_conv * (s_sink * s_sink)/((s_sink*s_sink - s_sp*s_sp) / 4);

	/*thermal conductances of thermal interface material*/
	thermo->pk_thermal->tim_gx = 1/getr(k_tim, cw, t_tim*ch);
	thermo->pk_thermal->tim_gy = 1/getr(k_tim, ch, t_tim*cw);
	thermo->pk_thermal->tim_gz = 1/getr(k_tim, t_tim, cw*ch);

}

/*fill in channel and wall conductances*/
void populate_channel_conductances(stk_descr_t *stk, thermal_data_t *thermo)
{
	fill_default_flowrate(stk, thermo);
	populate_h_ch(stk, thermo);

	int i,j,k;
	double cell_length=stk->length/stk->rows;
	int n_ch_per_layer=stk->microchannel->n_channels_per_layer;
	int n_ch_per_cell= n_ch_per_layer/stk->cols;
	double ch_width=stk->microchannel->channel_width;
	double ch_t=stk->microchannel->channel_height;
	double wall_width=stk->microchannel->wall_width;
	int nr=stk->rows;
	int chc_index; //channel cell index

	for(i=0;i<stk->n_layers-1;i++)
		for(j=0;j<nr;j++)
			for(k=0;k<n_ch_per_layer;k++)
		{
			chc_index=i*n_ch_per_layer*nr+j*n_ch_per_layer+k;
			/*liquid cell to its boudary conductances*/
			thermo->ch_thermal->ch_gx[chc_index]=ch_t*cell_length*thermo->ch_thermal->h_ch[chc_index];
			thermo->ch_thermal->ch_gz[chc_index]=ch_width*cell_length*thermo->ch_thermal->h_ch[chc_index];
		}

	/*wall cell to its x-adjacent wall cell conductances, only exist at boarder of each circuit cell.*/
	thermo->ch_thermal->wall_gx=2*K_SILICON*cell_length*ch_t/wall_width;
	/*wall cell to its y-adjacent wall cell conductances*/
	thermo->ch_thermal->wall_gy=K_SILICON*ch_t*wall_width/cell_length;
	/*wall cell to its z-boundary conductance*/
	thermo->ch_thermal->wall_gz=2*K_SILICON*cell_length*wall_width/ch_t;
}


/*calculate thermal transfer coefficient of each channel subject to their respective flowrate.*
 *Nusselt number: for entry length region, we use Seider-Tate correlation without consideration of viscosity change; otherwise Peng's correlation is used*
 *In practical, larger one is chosen*/
void populate_h_ch(stk_descr_t *stk, thermal_data_t *thermo)
{
	int i,j,k;
	double Nu_ST, Nu_P, Nu_Gz, *Re, Pr, Gz, DH,P, ch, alpha, beta, L0;
	int nl=stk->n_layers;
	int nr=stk->rows;
	int n_ch_per_layer=stk->microchannel->n_channels_per_layer;

	Re=(double *)malloc((nl-1)*n_ch_per_layer*sizeof(double));
	Pr=SP_FLUID*Dyn_Vis/(Density_Fluid*K_Fluid); //Prandtl Number
	P=2*(stk->microchannel->channel_height+stk->microchannel->channel_width); //Microchannel Perimeter
	DH=4*stk->microchannel->channel_height*stk->microchannel->channel_width/P; //Hydraulic diameter of microchannel
	ch=stk->length/nr; //cell height
	alpha=stk->microchannel->channel_height/stk->microchannel->channel_width; //channel aspect ratio
	beta=DH/(stk->microchannel->channel_width+stk->microchannel->wall_width); //channel density
	L0=0.015;

	double correction=4.5116;

	for(i=0;i<(nl-1)*n_ch_per_layer;i++)
		Re[i]=4*thermo->ch_thermal->flowrate[i]*Density_Fluid/(Dyn_Vis*P);

	printf("Re is %f, Pr is %f\n",Re[0], Pr);

	for(i=0;i<nl-1;i++)
		for(j=0;j<nr;j++)
			for(k=0;k<n_ch_per_layer;k++)
		{
			Nu_ST=0.1165*pow(beta, 0.81)*pow(alpha, -0.79)*pow(Re[i*n_ch_per_layer+k], 0.62)*pow(Pr, 0.333)
					+30*pow(Re[i*n_ch_per_layer+k], 0.333)*pow(Pr, 0.333)*pow(DH/((j+0.5)*ch), 1.67);

		//	Nu_ST=pow((pow((0.587*pow(((j+0.5)*ch/(DH*Re[i*n_ch_per_layer+k]*Pr)), -0.5)), 3)+558.4584), 1/3);
			thermo->ch_thermal->h_ch[i*n_ch_per_layer*nr+j*n_ch_per_layer+k]=K_Fluid*Nu_ST/DH;
		}

	free(Re);
}

void update_thermal_data(stk_descr_t *stk, thermal_data_t *thermo)
{
	populate_channel_conductances(stk, thermo);
}

void solve_steady_state(stk_descr_t *stk, thermal_data_t *thermo)
{
	if(stk->cooling_method==1)
		solve_steady_state_heatsink(stk, thermo);
	else if(stk->cooling_method==2)
		solve_steady_state_ch(stk, thermo);
	else
		fatal("undefined cooling method when trying to solve steady state temperatures.");
}

void solve_steady_state_heatsink(stk_descr_t *stk, thermal_data_t *thermo)
{
	int nl = stk->n_layers;
	int nr = stk->rows;
	int nc = stk->cols;

	int n_circuit=nl*nr*nc*4;
	int n_sp_core=nr*nc, n_sp_per=4;
	int n_sink_core=nr*nc, n_sink_per_inner=4, n_sink_per_outer=4;
	int size=n_circuit+n_sp_core+n_sp_per+n_sink_core+n_sink_per_inner+n_sink_per_outer;

	printf("n_circuit = %d, n_sp_core=%d, n_sink_core=%d\n", n_circuit, n_sp_core, n_sink_core);

	int *Ap, *Ai;
	double *Ax, *p, sum;
	int i,j;

	int layer, row, col;

	int *inc;

	Ap=(int *)calloc(size+1, sizeof(int));

	/*count for silicon circuit layers*/
	for(i=0;i<n_circuit;i++)
	{
		layer=double2int(i/(nr*nc));
		row=double2int(i/nc)%nr;
		col=i%nc;

		/*regarding the silicon circuit part, the matrix has non-zero elements for itself and adjacent silicon cells, take note of boundary cells.*/
		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1)  //south cell
			Ap[i+nc+1]++;
		if(layer!=0)  // cell above
			Ap[i-nc*nr+1]++;
		Ap[i+nc*nr+1]++; //cell below, spreader cell for layer closest to heatsink
	}

	/*count for heat spreader layer*/
	for(i=n_circuit;i<n_circuit+n_sp_core;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;

		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1) //south cell
			Ap[i+nc+1]++;

		// closest circuit layer
		Ap[i-nc*nr+1]++;
		//heat sink layer
		Ap[i+nc*nr+1]++;

		/*regarding the heat spreader periphery*/
		if(col==0) //west boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_W+1]++;
		if(col==nc-1) //east boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_E+1]++;
		if(row==0) //north boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_N+1]++;
		if(row==nr-1) //south boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_S+1]++;
	}

	/*heat sink layer*/
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;

		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1) //south cell
			Ap[i+nc+1]++;

		//top spreader layer
		Ap[i-nc*nr+1]++;

		/*regarding the heat sink periphery*/
		if(col==0) //west boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W+1]++;
		if(col==nc-1) //east boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E+1]++;
		if(row==0) //north boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N+1]++;
		if(row==nr-1) //south boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S+1]++;
	}

	/* heat spreader periphery */
	/* SP_W west spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_W+1]++; //diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W+1]++; //inner heat sink
	//heat spreader grid west boundary
	for(i=n_circuit;i<n_circuit+n_sp_core;i+=nc)
		Ap[i+1]++;
	/* SP_E east spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_E+1]++;//diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E+1]++; //inner heat sink
	//heat spreader grid east boundary
	for(i=n_circuit+nc-1;i<n_circuit+n_sp_core;i+=nc)
		Ap[i+1]++;
	/* SP_N north spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_N+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N+1]++;
	for(i=n_circuit;i<n_circuit+nc;i++)
		Ap[i+1]++;
	/*SP_S south spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_S+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S+1]++;
	for(i=n_circuit+n_sp_core-nc;i<n_circuit+n_sp_core;i++)
		Ap[i+1]++;

	/* heat sink periphery */
	/* SINK_C_W west inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W+1]++;//diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SP_W+1]++; //spreader
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_W+1]++; //outer heat sink
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i+=nc)
		Ap[i+1]++;
	/* SINK_C_E east inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E+1]++; //diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SP_E+1]++; //spreader
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_E+1]++; //outer heat sink
	for(i=n_circuit+n_sp_core+nc-1;i<n_circuit+n_sp_core+n_sink_core;i+=nc)
		Ap[i+1]++;
	/* SINK_C_N north inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N+1]++; //diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SP_N+1]++; //spreader
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_N+1]++;
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+nc;i++)
		Ap[i+1]++;
	/* SINK_C_S south inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SP_S+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_S+1]++;
	for(i=n_circuit+n_sp_core+n_sink_core-nc;i<n_circuit+n_sp_core+n_sink_core;i++)
		Ap[i+1]++;

	/* SINK_W west outer heat sink */
	for(i=n_circuit+n_sp_core+n_sink_core+SINK_W;i<size;i++)
	{
		Ap[i+1]++;
		Ap[i-3]++;
	}

	for(i=2;i<size+1;i++)
		Ap[i]+=Ap[i-1];

	printf("*************Ap Ready************\n");

	/*allocate space for Ai and Ax vector*/
	Ai=(int *)calloc(Ap[size], sizeof(int));
	Ax=(double *)calloc(Ap[size], sizeof(double));

	double *tsv_circuitx=stk->tsv_circuitx;
	double *tsv_circuity=stk->tsv_circuity;
	double *tsv_circuitz=stk->tsv_circuitz;
	double *tsv_bondx=stk->tsv_bondx;
	double *tsv_bondy=stk->tsv_bondy;
	double *tsv_bondz=stk->tsv_bondz;
	package_thermal_t *pk=thermo->pk_thermal;
	double rz_BEOL, rz_source, rz_bulk, rz_TIM;
	double k_tsv=stk->k_tsv, k_TIM=stk->heatsink->k_TIM;

	inc=(int *)calloc(size, sizeof(int));

	for(i=0;i<size;i++)
		inc[i]=0;

	int diag_idx;



	/*prepare Ai and Ax vector*/
	//prepare fine-grained layer thermal parameters
	double *layer_gx=(double *)malloc(nl*4*sizeof(double));
	double *layer_gy=(double *)malloc(nl*4*sizeof(double));
	double *layer_gz=(double *)malloc(nl*4*sizeof(double));
	double *layer_k=(double *)malloc(nl*4*sizeof(double));
	double *tsvx=(double *)malloc((nl*4-1)*nr*nc*sizeof(double));
	double *tsvy=(double *)malloc((nl*4-1)*nr*nc*sizeof(double));
	double *tsvz=(double *)malloc((nl*4-1)*nr*nc*sizeof(double));
	for(i=0;i<nl;i++)
	{
		layer_gx[i*4]=thermo->gx_BEOL[i];
		layer_gx[i*4+1]=thermo->gx_source[i];
		layer_gx[i*4+2]=thermo->gx_bulk[i];
		layer_gx[i*4+3]=pk->tim_gx;
	}
	for(i=0;i<nl;i++)
	{
		layer_gy[i*4]=thermo->gy_BEOL[i];
		layer_gy[i*4+1]=thermo->gy_source[i];
		layer_gy[i*4+2]=thermo->gy_bulk[i];
		layer_gy[i*4+3]=pk->tim_gy;
	}
	for(i=0;i<nl;i++)
	{
		layer_gz[i*4]=thermo->gz_BEOL[i];
		layer_gz[i*4+1]=thermo->gz_source[i];
		layer_gz[i*4+2]=thermo->gz_bulk[i];
		layer_gz[i*4+3]=pk->tim_gz;
	}
	for(i=0;i<nl;i++)
	{
		layer_k[i*4]=K_BEOL;
		layer_k[i*4+1]=K_SILICON;
		layer_k[i*4+2]=K_SILICON;
		layer_k[i*4+3]=k_TIM;
	}
	for(i=0;i<nl;i++)
	{
		for(j=0;j<3*nr*nc;j++)
		{
			tsvx[i*4*nr*nc+j]=stk->tsv_circuitx[i*3*nr*nc+j];
			tsvy[i*4*nr*nc+j]=stk->tsv_circuity[i*3*nr*nc+j];
			tsvz[i*4*nr*nc+j]=stk->tsv_circuitz[i*3*nr*nc+j];
		}
		if(i!=nl-1)
			for(j=0;j<nr*nc;j++)
			{
				tsvx[i*4*nr*nc+3*nr*nc+j]=stk->tsv_bondx[i*nr*nc+j];
				tsvy[i*4*nr*nc+3*nr*nc+j]=stk->tsv_bondy[i*nr*nc+j];
				tsvz[i*4*nr*nc+3*nr*nc+j]=stk->tsv_bondz[i*nr*nc+j];
			}
	}

	//circuit layers first
	for(i=0;i<n_circuit;i++)
	{
		layer=double2int(i/(nr*nc));
		row=double2int(i/nc)%nr;
		col=i%nc;

		sum=0;
		/*regarding the silicon circuit part, the matrix has non-zero elements for itself and adjacent silicon cells, take note of boundary cells.*/
		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //conductance to west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-(1+0.5*(tsvx[i]+tsvx[i-1]))*layer_gx[layer];
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}

		if(col!=nc-1) //conductance to east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-(1+0.5*(tsvx[i]+tsvx[i+1]))*layer_gx[layer];
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}

		if(row!=0) //conductance to north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-(1+0.5*(tsvy[i]+tsvy[i-nc]))*layer_gy[layer];
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}

		if(row!=nr-1) //conductance to south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-(1+0.5*(tsvy[i]+tsvy[i+nc]))*layer_gy[layer];
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}

		if(layer!=0) //conductance to top layer except top most layer
		{
			Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
			Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-(1+tsvz[i-nc*nr])*layer_gz[layer-1];
			sum-=Ax[Ap[i-nc*nr]+(inc[i-nc*nr]-1)];
		}

		//conductance to bottom layer
		Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
		if(layer!=nl*4-1)
			Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-(1+tsvz[i])*layer_gz[layer];
		else
			Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-layer_gz[layer];
		sum-=Ax[Ap[i+nc*nr]+(inc[i+nc*nr]-1)];


		Ax[diag_idx]=sum;

	}


	/* heat spreader layer core*/
	for(i=n_circuit;i<n_circuit+n_sp_core;i++)  //heat spreader core part
	{
		row=double2int(i/nc)%nr;
		col=i%nc;

		sum=0;

		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //conductance to west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-1/pk->sp_rx;
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++)]=
					-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]-1];
		}

		if(col!=nc-1) //conductance to east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-1/pk->sp_rx;
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++)]=
					-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]-1];
		}

		if(row!=0) //conductance to north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-1/pk->sp_ry;
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++)]=
					-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]-1];
		}

		if(row!=nr-1) //conductance to south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-1/pk->sp_ry;
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++)]=
					-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]-1];
		}

		// circuit layer next to heat spreader
		Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
/*		rz_source=0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[((nl-1)*NumSubLayers+1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[nl-1]);
		rz_bulk=1.0/((1+(k_tsv-K_SILICON)*tsv_circuit[((nl-1)*NumSubLayers+2)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_bulk[nl-1]);
		rz_TIM=1.0/pk->tim_gz;
		Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(rz_source+rz_bulk+rz_TIM+pk->sp_rz);*/

		Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-layer_gz[nl*4-1];
		sum-=Ax[Ap[i-nc*nr]+inc[i-nc*nr]-1];

		// heat sink layer cell next to heat spreader
		Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
		//Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(pk->sp_rz+pk->sink_rz);
		Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(2*pk->sp_rz);
		sum-=Ax[Ap[i+nc*nr]+inc[i+nc*nr]-1];

		Ax[diag_idx]=sum;
	}

	/* heat sink layer */
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;
		sum=0;

		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //conductance to west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-1/pk->sink_rx;
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=
					-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]-1];
		}

		if(col!=nc-1) //conductance to east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-1/pk->sink_rx;
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=
					-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]-1];
		}

		if(row!=0) //conductance to north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-1/pk->sink_ry;
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=
					-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]-1];
		}

		if(row!=nr-1) //conductance to south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-1/pk->sink_ry;
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=
					-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]-1];
		}

		// heat spreader layer cell
		Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
		//Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(pk->sp_rz+pk->sink_rz);
		Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(2*pk->sp_rz);
		sum-=Ax[Ap[i-nc*nr]+inc[i-nc*nr]-1];

		// heat sink is connected to ambient
		//sum+=1.0/(pk->sink_rz+pk->sink_ramb);
		sum+=1.0/(2*pk->sink_rz+pk->sink_ramb);
		Ax[diag_idx]=sum;
	}

	/* SP_W west spreader periphery */

	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_W;
	Ax[diag_idx]=0;

	/*to spreader core*/
	for(i=n_circuit;i<n_circuit+n_sp_core;i+=nc)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_W;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]]=n_circuit+n_sp_core+n_sink_core+SP_W;
	//Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]-1];

	Ax[diag_idx]=sum;

	/* SP_E east spreader periphery */
	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_E;
	Ax[diag_idx]=0;

	for(i=n_circuit+nc-1;i<n_circuit+n_sp_core;i+=nc)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_E;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]]=n_circuit+n_sp_core+n_sink_core+SP_E;
	//Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]-1];


	Ax[diag_idx]=sum;

	/* SP_N north spreader periphery */
	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_N;
	Ax[diag_idx]=0;

	/*to spreader core*/
	for(i=n_circuit;i<n_circuit+nc;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_N;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]]=n_circuit+n_sp_core+n_sink_core+SP_N;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]-1];

	Ax[diag_idx]=sum;

	/* SP_S south spreader periphery */
	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_S;
	Ax[diag_idx]=0;

	/*to spreader core*/
	for(i=n_circuit+n_sp_core-nc;i<n_circuit+n_sp_core;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_S;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]]=n_circuit+n_sp_core+n_sink_core+SP_S;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]-1];

	Ax[diag_idx]=sum;

	/* SINK_C_W west inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
	Ax[diag_idx]=0;

	//to spreader peripheral
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i+=nc)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_W]]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_W]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_W]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_we+pk->ramb_c_we);
	sum+=1.0/(2*pk->rz_hs_c_we+pk->ramb_c_we);
	Ax[diag_idx]=sum;

	/* SINK_C_E east inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
	Ax[diag_idx]=0;

	//spreader
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core+nc-1;i<n_circuit+n_sp_core+n_sink_core;i+=nc) {
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_E]]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_E]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_E]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_we+pk->ramb_c_we);
	sum+=1.0/(2*pk->rz_hs_c_we+pk->ramb_c_we);
	Ax[diag_idx]=sum;

	/* SINK_C_N north inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
	Ax[diag_idx]=0;

	//spreader
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+nc;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_N]]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_N]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_N]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_ns+pk->ramb_c_ns);
	sum+=1.0/(2*pk->rz_hs_c_ns+pk->ramb_c_ns);
	Ax[diag_idx]=sum;

	/* SINK_C_S south inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
	Ax[diag_idx]=0;

	//spreader
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core+n_sink_core-nc;i<n_circuit+n_sp_core+n_sink_core;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_S]]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_S]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_S]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_ns+pk->ramb_c_ns);
	sum+=1.0/(2*pk->rz_hs_c_ns+pk->ramb_c_ns);
	Ax[diag_idx]=sum;

	/* SINK_W west outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_W]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_W;
	Ax[diag_idx]=0;

	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]]=n_circuit+n_sp_core+n_sink_core+SINK_W;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	/* SINK_E east outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_E]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_E;
	Ax[diag_idx]=0;
	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]]=n_circuit+n_sp_core+n_sink_core+SINK_E;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]-1];
	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	/* SINK_N north outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_N]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_N;
	Ax[diag_idx]=0;
	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]]=n_circuit+n_sp_core+n_sink_core+SINK_N;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]-1];
	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	/* SINK_S south outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_S]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_S;
	Ax[diag_idx]=0;
	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]]=n_circuit+n_sp_core+n_sink_core+SINK_S;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]-1];
	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	printf("*************Ai and Ax Ready************\n");

	/*initialize power matrix*/

	p=(double *)calloc(size, sizeof(double));

	for(i=0;i<size;i++)
		p[i]=0;

	for(i=0;i<nl;i++)
		for(j=0;j<nr*nc;j++)
			p[(i*4+1)*nr*nc+j]=thermo->power_grid[i*nr*nc+j];

	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i++)
		p[i]+=AMBIENT/(2*pk->sink_rz+pk->sink_ramb);

	for(i=n_circuit+n_sp_core+n_sink_core+SINK_C_W;i<n_circuit+n_sp_core+n_sink_core+SINK_C_E+1;i++)
		p[i]+=AMBIENT/(2*pk->rz_hs_c_we+pk->ramb_c_we);

	for(i=n_circuit+n_sp_core+n_sink_core+SINK_C_N;i<n_circuit+n_sp_core+n_sink_core+SINK_C_S+1;i++)
		p[i]+=AMBIENT/(2*pk->rz_hs_c_ns+pk->ramb_c_ns);

	for(i=n_circuit+n_sp_core+n_sink_core+SINK_W;i<size;i++)
		p[i]+=AMBIENT/(2*pk->rz_hs_outer+pk->ramb_outer);

/*	for(i=0;i<size;i++)
		printf("%f\n", p[i]);
*/
	printf("*************p Ready************\n");

	//	print_klu_matrix(Ap, Ai, Ax, p, size);

	klu_symbolic *Symbolic ;
	klu_numeric *Numeric ;
	klu_common Common ;

	// printf("*************KLU variable declared************\n");
	klu_defaults (&Common) ;
	Symbolic = klu_analyze (size, Ap, Ai, &Common) ;
	//printf("*************Symbolic Ready************\n");
	Numeric = klu_factor (Ap, Ai, Ax, Symbolic, &Common) ;
	//printf("*************Numeric Ready************\n");
	klu_solve (Symbolic, Numeric, size, 1, p, &Common) ;
	printf("*************KLU solved************\n");
	klu_free_symbolic (&Symbolic, &Common) ;
	klu_free_numeric (&Numeric, &Common) ;
	/*output matrix elements to a file, for debug purpose*/
	//  dump_sparse_matrix(Ap, Ai, Ax, size, "matrixdetail");


/*	for(i=0;i<size;++i)
	{
		printf("%d :\t", i);
		for(j=Ap[i];j<Ap[i+1];j++)
			printf("%lf\t", Ax[j]);
		printf("\n");
	}*/
	free(Ai);
	free(Ap);
	free(Ax);

	for(i=0;i<nl;i++)
		for(j=0;j<nr*nc;j++)
		thermo->steady_temp_grid[i*nr*nc+j]=p[(i*4+1)*nr*nc+j];



	/*map back to n units*/
	xlate_temp_g2b(stk,thermo);
	free(p);

}

/*compute steady state temperature under heatsink cooling*/
void solve_steady_state_heatsink_Coarse(stk_descr_t *stk, thermal_data_t *thermo)
{
	int nl = stk->n_layers;
	int nr = stk->rows;
	int nc = stk->cols;

	int n_circuit=nl*nr*nc;
	int n_sp_core=nr*nc, n_sp_per=4;
	int n_sink_core=nr*nc, n_sink_per_inner=4, n_sink_per_outer=4;
	int size=n_circuit+n_sp_core+n_sp_per+n_sink_core+n_sink_per_inner+n_sink_per_outer;

	printf("n_circuit = %d, n_sp_core=%d, n_sink_core=%d\n", n_circuit, n_sp_core, n_sink_core);

	int *Ap, *Ai;
	double *Ax, *p, sum;
	int i,j;

	int layer, row, col;

	int *inc;

	Ap=(int *)calloc(size+1, sizeof(int));

	/*count for silicon circuit layers*/
	for(i=0;i<n_circuit;i++)
	{
		layer=double2int(i/(nr*nc));
		row=double2int(i/nc)%nr;
		col=i%nc;

		/*regarding the silicon circuit part, the matrix has non-zero elements for itself and adjacent silicon cells, take note of boundary cells.*/
		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1)  //south cell
			Ap[i+nc+1]++;
		if(layer!=0)  // cell above
			Ap[i-nc*nr+1]++;
		Ap[i+nc*nr+1]++; //cell below, spreader cell for layer closest to heatsink
	}

	/*count for heat spreader layer*/
	for(i=n_circuit;i<n_circuit+n_sp_core;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;

		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1) //south cell
			Ap[i+nc+1]++;

		// closest circuit layer
		Ap[i-nc*nr+1]++;
		//heat sink layer
		Ap[i+nc*nr+1]++;

		/*regarding the heat spreader periphery*/
		if(col==0) //west boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_W+1]++;
		if(col==nc-1) //east boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_E+1]++;
		if(row==0) //north boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_N+1]++;
		if(row==nr-1) //south boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SP_S+1]++;
	}

	/*heat sink layer*/
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;

		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1) //south cell
			Ap[i+nc+1]++;

		//top spreader layer
		Ap[i-nc*nr+1]++;

		/*regarding the heat sink periphery*/
		if(col==0) //west boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W+1]++;
		if(col==nc-1) //east boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E+1]++;
		if(row==0) //north boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N+1]++;
		if(row==nr-1) //south boundary
			Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S+1]++;
	}

	/* heat spreader periphery */
	/* SP_W west spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_W+1]++; //diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W+1]++; //inner heat sink
	//heat spreader grid west boundary
	for(i=n_circuit;i<n_circuit+n_sp_core;i+=nc)
		Ap[i+1]++;
	/* SP_E east spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_E+1]++;//diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E+1]++; //inner heat sink
	//heat spreader grid east boundary
	for(i=n_circuit+nc-1;i<n_circuit+n_sp_core;i+=nc)
		Ap[i+1]++;
	/* SP_N north spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_N+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N+1]++;
	for(i=n_circuit;i<n_circuit+nc;i++)
		Ap[i+1]++;
	/*SP_S south spreader periphery */
	Ap[n_circuit+n_sp_core+n_sink_core+SP_S+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S+1]++;
	for(i=n_circuit+n_sp_core-nc;i<n_circuit+n_sp_core;i++)
		Ap[i+1]++;

	/* heat sink periphery */
	/* SINK_C_W west inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W+1]++;//diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SP_W+1]++; //spreader
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_W+1]++; //outer heat sink
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i+=nc)
		Ap[i+1]++;
	/* SINK_C_E east inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E+1]++; //diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SP_E+1]++; //spreader
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_E+1]++; //outer heat sink
	for(i=n_circuit+n_sp_core+nc-1;i<n_circuit+n_sp_core+n_sink_core;i+=nc)
		Ap[i+1]++;
	/* SINK_C_N north inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N+1]++; //diagonal
	Ap[n_circuit+n_sp_core+n_sink_core+SP_N+1]++; //spreader
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_N+1]++;
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+nc;i++)
		Ap[i+1]++;
	/* SINK_C_S south inner heat sink */
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SP_S+1]++;
	Ap[n_circuit+n_sp_core+n_sink_core+SINK_S+1]++;
	for(i=n_circuit+n_sp_core+n_sink_core-nc;i<n_circuit+n_sp_core+n_sink_core;i++)
		Ap[i+1]++;

	/* SINK_W west outer heat sink */
	for(i=n_circuit+n_sp_core+n_sink_core+SINK_W;i<size;i++)
	{
		Ap[i+1]++;
		Ap[i-3]++;
	}

	for(i=2;i<size+1;i++)
		Ap[i]+=Ap[i-1];

	printf("*************Ap Ready************\n");

	/*allocate space for Ai and Ax vector*/
	Ai=(int *)calloc(Ap[size], sizeof(int));
	Ax=(double *)calloc(Ap[size], sizeof(double));

	double *tsv_circuitx=stk->tsv_circuitx;
	double *tsv_circuity=stk->tsv_circuity;
	double *tsv_circuitz=stk->tsv_circuitz;
	double *tsv_bondx=stk->tsv_bondx;
	double *tsv_bondy=stk->tsv_bondy;
	double *tsv_bondz=stk->tsv_bondz;

	package_thermal_t *pk=thermo->pk_thermal;
	double rz_BEOL, rz_source, rz_bulk, rz_TIM;
	double k_tsv=stk->k_tsv, k_TIM=stk->heatsink->k_TIM;

	inc=(int *)calloc(size, sizeof(int));

	for(i=0;i<size;i++)
		inc[i]=0;

	int diag_idx;



	/*prepare Ai and Ax vector*/

	//circuit layers first
	for(i=0;i<n_circuit;i++)
	{
		layer=double2int(i/(nr*nc));
		row=double2int(i/nc)%nr;
		col=i%nc;

		sum=0;
		/*regarding the silicon circuit part, the matrix has non-zero elements for itself and adjacent silicon cells, take note of boundary cells.*/
		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //conductance to west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-(thermo->gx_BEOL[layer]+thermo->gx_source[layer]+thermo->gx_bulk[layer]+pk->tim_gx);
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}

		if(col!=nc-1) //conductance to east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-(thermo->gx_BEOL[layer]+thermo->gx_source[layer]+thermo->gx_bulk[layer]+pk->tim_gx);
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}

		if(row!=0) //conductance to north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-(thermo->gy_BEOL[layer]+thermo->gy_source[layer]+thermo->gy_bulk[layer]+pk->tim_gy);
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}

		if(row!=nr-1) //conductance to south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-(thermo->gy_BEOL[layer]+thermo->gy_source[layer]+thermo->gy_bulk[layer]+pk->tim_gy);
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}

/*		if(layer!=0) //conductance to top layer except top most layer
		{
			Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
			rz_BEOL=1.0/((1+(k_tsv-K_BEOL)*tsv_circuit[layer*NumSubLayers*nr*nc+row*nc+col]/K_BEOL)*thermo->gz_BEOL[layer]);
			rz_source=0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers+1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[layer])
					+0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers-2)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[layer-1]);
			rz_bulk=1.0/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers-1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_bulk[layer-1]);
			rz_TIM=1.0/((1+(k_tsv-k_TIM)*tsv_bond[(layer-1)*nr*nc+row*nc+col]/k_TIM)*pk->tim_gz);

			Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(rz_BEOL+rz_source+rz_bulk+rz_TIM);
			sum-=Ax[Ap[i-nc*nr]+(inc[i-nc*nr]-1)];
		}
*/
		if(layer!=0) //conductance to top layer except top most layer
		{
			Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
			rz_BEOL=1.0/thermo->gz_BEOL[layer-1];
			rz_source=1/thermo->gz_source[layer-1];
			rz_bulk=1.0/thermo->gz_bulk[layer-1];
			rz_TIM=1.0/pk->tim_gz;

			Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(rz_BEOL+rz_source+rz_bulk+rz_TIM);
			sum-=Ax[Ap[i-nc*nr]+(inc[i-nc*nr]-1)];
		}



/*		if(layer!=nl-1) //conductance to bottom layer except the bottom most layer
		{
			Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
			rz_BEOL=1.0/((1+(k_tsv-K_BEOL)*tsv_circuit[(layer+1)*NumSubLayers*nr*nc+row*nc+col]/K_BEOL)*thermo->gz_BEOL[layer+1]);
			rz_source=0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers+1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[layer])
					+0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[((layer+1)*NumSubLayers+1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[layer+1]);
			rz_bulk=1.0/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers+2)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_bulk[layer]);
			rz_TIM=1.0/((1+(k_tsv-k_TIM)*tsv_bond[layer*nr*nc+row*nc+col]/k_TIM)*pk->tim_gz);

			Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(rz_BEOL+rz_source+rz_bulk+rz_TIM);
			sum-=Ax[Ap[i+nc*nr]+(inc[i+nc*nr]-1)];
		}
		else if(layer==nl-1) // bottom most layer next to heat spreader
		{
			Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
			rz_source=0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers+1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[layer]);
			rz_bulk=1.0/((1+(k_tsv-K_SILICON)*tsv_circuit[(layer*NumSubLayers+2)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_bulk[layer]);
			rz_TIM=1.0/thermo->pk_thermal->tim_gz;

			Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(rz_source+rz_bulk+rz_TIM+pk->sp_rz);
			sum-=Ax[Ap[i+nc*nr]+(inc[i+nc*nr]-1)];
		}
*/

//		if(layer!=nl-1)
//		{
			Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
			rz_BEOL=1.0/thermo->gz_BEOL[layer];
			rz_source=1/thermo->gz_source[layer];
			rz_bulk=1.0/thermo->gz_bulk[layer];
			rz_TIM=1.0/pk->tim_gz;
			Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(rz_BEOL+rz_source+rz_bulk+rz_TIM);
			sum-=Ax[Ap[i+nc*nr]+(inc[i+nc*nr]-1)];
//		}

	/*	else
		{
			Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
			Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(2*pk->sp_rz);
			sum-=Ax[Ap[i+nc*nr]+(inc[i+nc*nr]-1)];
		}
*/

		Ax[diag_idx]=sum;

	}


	/* heat spreader layer core*/
	for(i=n_circuit;i<n_circuit+n_sp_core;i++)  //heat spreader core part
	{
		row=double2int(i/nc)%nr;
		col=i%nc;

		sum=0;

		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //conductance to west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-1/pk->sp_rx;
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++)]=
					-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]-1];
		}

		if(col!=nc-1) //conductance to east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-1/pk->sp_rx;
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++)]=
					-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]-1];
		}

		if(row!=0) //conductance to north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-1/pk->sp_ry;
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++)]=
					-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]-1];
		}

		if(row!=nr-1) //conductance to south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-1/pk->sp_ry;
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++)]=
					-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]-1];
		}

		// circuit layer next to heat spreader
		Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
/*		rz_source=0.5/((1+(k_tsv-K_SILICON)*tsv_circuit[((nl-1)*NumSubLayers+1)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_source[nl-1]);
		rz_bulk=1.0/((1+(k_tsv-K_SILICON)*tsv_circuit[((nl-1)*NumSubLayers+2)*nr*nc+row*nc+col]/K_SILICON)*thermo->gz_bulk[nl-1]);
		rz_TIM=1.0/pk->tim_gz;
		Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(rz_source+rz_bulk+rz_TIM+pk->sp_rz);*/

		rz_BEOL=1/thermo->gz_BEOL[nl-1];
		rz_source=1/thermo->gz_source[nl-1];
		rz_bulk=1/thermo->gz_bulk[nl-1];
		rz_TIM=1/pk->tim_gz;
		Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(rz_BEOL+rz_source+rz_bulk+rz_TIM);
		sum-=Ax[Ap[i-nc*nr]+inc[i-nc*nr]-1];

		// heat sink layer cell next to heat spreader
		Ai[Ap[i+nc*nr]+inc[i+nc*nr]]=i;
		//Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(pk->sp_rz+pk->sink_rz);
		Ax[Ap[i+nc*nr]+(inc[i+nc*nr]++)]=-1.0/(2*pk->sp_rz);
		sum-=Ax[Ap[i+nc*nr]+inc[i+nc*nr]-1];

		Ax[diag_idx]=sum;
	}

	/* heat sink layer */
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;
		sum=0;

		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //conductance to west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-1/pk->sink_rx;
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=
					-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]-1];
		}

		if(col!=nc-1) //conductance to east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-1/pk->sink_rx;
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=
					-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]-1];
		}

		if(row!=0) //conductance to north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-1/pk->sink_ry;
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=
					-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]-1];
		}

		if(row!=nr-1) //conductance to south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-1/pk->sink_ry;
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}
		else
		{
			Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]]=i;
			Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=
					-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
			sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]-1];
		}

		// heat spreader layer cell
		Ai[Ap[i-nc*nr]+inc[i-nc*nr]]=i;
		//Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(pk->sp_rz+pk->sink_rz);
		Ax[Ap[i-nc*nr]+(inc[i-nc*nr]++)]=-1.0/(2*pk->sp_rz);
		sum-=Ax[Ap[i-nc*nr]+inc[i-nc*nr]-1];

		// heat sink is connected to ambient
		//sum+=1.0/(pk->sink_rz+pk->sink_ramb);
		sum+=1.0/(2*pk->sink_rz+pk->sink_ramb);
		Ax[diag_idx]=sum;
	}

	/* SP_W west spreader periphery */

	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_W;
	Ax[diag_idx]=0;

	/*to spreader core*/
	for(i=n_circuit;i<n_circuit+n_sp_core;i+=nc)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_W;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]]=n_circuit+n_sp_core+n_sink_core+SP_W;
	//Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]-1];

	Ax[diag_idx]=sum;

	/* SP_E east spreader periphery */
	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_E;
	Ax[diag_idx]=0;

	for(i=n_circuit+nc-1;i<n_circuit+n_sp_core;i+=nc)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_E;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_rx+nr*pk->rx_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]]=n_circuit+n_sp_core+n_sink_core+SP_E;
	//Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]-1];


	Ax[diag_idx]=sum;

	/* SP_N north spreader periphery */
	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_N;
	Ax[diag_idx]=0;

	/*to spreader core*/
	for(i=n_circuit;i<n_circuit+nc;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_N;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]]=n_circuit+n_sp_core+n_sink_core+SP_N;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]-1];

	Ax[diag_idx]=sum;

	/* SP_S south spreader periphery */
	sum=0;
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SP_S;
	Ax[diag_idx]=0;

	/*to spreader core*/
	for(i=n_circuit+n_sp_core-nc;i<n_circuit+n_sp_core;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SP_S;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sp_ry+nc*pk->ry_sp_per);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	/*to heat sink peripheral*/
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]]=n_circuit+n_sp_core+n_sink_core+SP_S;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]-1];

	Ax[diag_idx]=sum;

	/* SINK_C_W west inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
	Ax[diag_idx]=0;

	//to spreader peripheral
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+(inc[n_circuit+n_sp_core+n_sink_core+SP_W]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_W]+inc[n_circuit+n_sp_core+n_sink_core+SP_W]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i+=nc)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_W]]=n_circuit+n_sp_core+n_sink_core+SINK_C_W;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_W]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_W]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_we+pk->ramb_c_we);
	sum+=1.0/(2*pk->rz_hs_c_we+pk->ramb_c_we);
	Ax[diag_idx]=sum;

	/* SINK_C_E east inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
	Ax[diag_idx]=0;

	//spreader
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++)]=-1/(pk->rz_sp_we+pk->rz_hs_c_we);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+(inc[n_circuit+n_sp_core+n_sink_core+SP_E]++)]=-1/(2*pk->rz_sp_we);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_E]+inc[n_circuit+n_sp_core+n_sink_core+SP_E]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core+nc-1;i<n_circuit+n_sp_core+n_sink_core;i+=nc) {
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_rx+nr*pk->rx_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_E]]=n_circuit+n_sp_core+n_sink_core+SINK_C_E;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_E]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_E]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_we+pk->ramb_c_we);
	sum+=1.0/(2*pk->rz_hs_c_we+pk->ramb_c_we);
	Ax[diag_idx]=sum;

	/* SINK_C_N north inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
	Ax[diag_idx]=0;

	//spreader
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+(inc[n_circuit+n_sp_core+n_sink_core+SP_N]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_N]+inc[n_circuit+n_sp_core+n_sink_core+SP_N]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+nc;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_N]]=n_circuit+n_sp_core+n_sink_core+SINK_C_N;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_N]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_N]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_ns+pk->ramb_c_ns);
	sum+=1.0/(2*pk->rz_hs_c_ns+pk->ramb_c_ns);
	Ax[diag_idx]=sum;

	/* SINK_C_S south inner heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
	Ax[diag_idx]=0;

	//spreader
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
//	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++)]=-1/(pk->rz_sp_ns+pk->rz_hs_c_ns);
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+(inc[n_circuit+n_sp_core+n_sink_core+SP_S]++)]=-1/(2*pk->rz_sp_ns);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SP_S]+inc[n_circuit+n_sp_core+n_sink_core+SP_S]-1];

	//heat sink core
	for(i=n_circuit+n_sp_core+n_sink_core-nc;i<n_circuit+n_sp_core+n_sink_core;i++)
	{
		Ai[Ap[i]+inc[i]]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
		Ax[Ap[i]+(inc[i]++)]=-1.0/(0.5*pk->sink_ry+nc*pk->ry_hs_ci);
		sum-=Ax[Ap[i]+inc[i]-1];
	}

	//outer heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_S]]=n_circuit+n_sp_core+n_sink_core+SINK_C_S;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_S]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_S]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_c_ns+pk->ramb_c_ns);
	sum+=1.0/(2*pk->rz_hs_c_ns+pk->ramb_c_ns);
	Ax[diag_idx]=sum;

	/* SINK_W west outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_W]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_W;
	Ax[diag_idx]=0;

	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]]=n_circuit+n_sp_core+n_sink_core+SINK_W;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_W]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_W]-1];

	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	/* SINK_E east outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_E]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_E;
	Ax[diag_idx]=0;
	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]]=n_circuit+n_sp_core+n_sink_core+SINK_E;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]++)]=-1.0/(pk->rx_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_E]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_E]-1];
	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	/* SINK_N north outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_N]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_N;
	Ax[diag_idx]=0;
	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]]=n_circuit+n_sp_core+n_sink_core+SINK_N;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_N]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_N]-1];
	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	/* SINK_S south outer heat sink */
	sum=0;
	//diagonal
	diag_idx=Ap[n_circuit+n_sp_core+n_sink_core+SINK_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_S]++);
	Ai[diag_idx]=n_circuit+n_sp_core+n_sink_core+SINK_S;
	Ax[diag_idx]=0;
	//inner heat sink
	Ai[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]]=n_circuit+n_sp_core+n_sink_core+SINK_S;
	Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+(inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]++)]=-1.0/(pk->ry_hs_co+pk->rxy_hs_outer);
	sum-=Ax[Ap[n_circuit+n_sp_core+n_sink_core+SINK_C_S]+inc[n_circuit+n_sp_core+n_sink_core+SINK_C_S]-1];
	//ambient
//	sum+=1.0/(pk->rz_hs_outer+pk->ramb_outer);
	sum+=1.0/(2*pk->rz_hs_outer+pk->ramb_outer);
	Ax[diag_idx]=sum;

	printf("*************Ai and Ax Ready************\n");

	/*initialize power matrix*/

	p=(double *)calloc(size, sizeof(double));

	for(i=0;i<n_circuit;i++)
		p[i]=thermo->power_grid[i];

	for(i=n_circuit;i<size;++i)
		p[i]=0;

	for(i=n_circuit+n_sp_core;i<n_circuit+n_sp_core+n_sink_core;i++)
		p[i]+=AMBIENT/(2*pk->sink_rz+pk->sink_ramb);

	for(i=n_circuit+n_sp_core+n_sink_core+SINK_C_W;i<n_circuit+n_sp_core+n_sink_core+SINK_C_E+1;i++)
		p[i]+=AMBIENT/(2*pk->rz_hs_c_we+pk->ramb_c_we);

	for(i=n_circuit+n_sp_core+n_sink_core+SINK_C_N;i<n_circuit+n_sp_core+n_sink_core+SINK_C_S+1;i++)
		p[i]+=AMBIENT/(2*pk->rz_hs_c_ns+pk->ramb_c_ns);

	for(i=n_circuit+n_sp_core+n_sink_core+SINK_W;i<size;i++)
		p[i]+=AMBIENT/(2*pk->rz_hs_outer+pk->ramb_outer);

	printf("*************p Ready************\n");

	//	print_klu_matrix(Ap, Ai, Ax, p, size);

	klu_symbolic *Symbolic ;
	klu_numeric *Numeric ;
	klu_common Common ;

	// printf("*************KLU variable declared************\n");
	klu_defaults (&Common) ;
	Symbolic = klu_analyze (size, Ap, Ai, &Common) ;
	//printf("*************Symbolic Ready************\n");
	Numeric = klu_factor (Ap, Ai, Ax, Symbolic, &Common) ;
	//printf("*************Numeric Ready************\n");
	klu_solve (Symbolic, Numeric, size, 1, p, &Common) ;
	printf("*************KLU solved************\n");
	klu_free_symbolic (&Symbolic, &Common) ;
	klu_free_numeric (&Numeric, &Common) ;
	/*output matrix elements to a file, for debug purpose*/
	//  dump_sparse_matrix(Ap, Ai, Ax, size, "matrixdetail");


	for(i=0;i<Ap[size];++i){
		printf("%lf ", Ax[i]);
		if(i%10==9)printf("\n");
	}
	free(Ai);
	free(Ap);
	free(Ax);

	for(i=0;i<n_circuit;i++)
	{
		thermo->steady_temp_grid[i]=p[i];
	}

	for(i=size-8;i<size;i++)
		printf("%f\n", p[i]);
	/*map back to n units*/
	xlate_temp_g2b(stk,thermo);
	free(p);

}

/*compute steady state temperature with microfluidic cooling*/
void solve_steady_state_ch(stk_descr_t *stk, thermal_data_t *thermo)
{

	int nl = stk->n_layers;
	int nr = stk->rows;
	int nc = stk->cols;
	int n_ch_per_layer=stk->microchannel->n_channels_per_layer;
	int n_ch_per_cell= n_ch_per_layer/nc;

	int i,j,*inc, ch_num, diag_idx, chc_index,index,base;
	int layer, row, col, layer1, row1, col1;
	int n_circuit=nl*nr*nc, nwall=(n_ch_per_cell+1)*nr*nc*(nl-1), nch=(nl-1)*n_ch_per_layer*nr;
	int size = n_circuit+nwall+nch;


	printf("matrix size is %d x %d, n_circuit=%d, nwall=%d, nch=%d\n", size, size, n_circuit, nwall, nch);

	int *Ap, *Ai;
	double *Ax, *p, sum;

	Ap=(int *)calloc(size+1, sizeof(int));

	/*prepare Ap vector, decide quantity of nonzeros row by row*/
	for(i=0;i<size+1;i++)
		Ap[i]=0;

	/*count for silicon circuit layers*/
	for(i=0;i<n_circuit;i++)
	{
		layer=double2int(i/(nr*nc));
		row=double2int(i/nc)%nr;
		col=i%nc;

		/*regarding the silicon circuit part, the matrix has non-zero elements for itself and adjacent silicon cells, take note of boundary cells.*/
		Ap[i+1]++; //diagonal

		if(col!=0) //west cell
			Ap[i]++;
		if(col!=nc-1) //east cell
			Ap[i+2]++;
		if(row!=0) //north cell
			Ap[i-nc+1]++;
		if(row!=nr-1) //south cell
			Ap[i+nc+1]++;

		/*regarding the channel wall part, the matrix has non-zero elements for wall cells above and below.*/
		if(layer!=nl-1) //wall cell below
			for(j=0;j<n_ch_per_cell+1;j++)
				Ap[n_circuit+i*(n_ch_per_cell+1)+j+1]++;

		if(layer!=0) //wall cell above
			for(j=0;j<n_ch_per_cell+1;j++)
				Ap[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j+1]++;

		/*regarding the channel part, the matrix has non-zero elements for channel cells above and below.*/
		if(layer!=nl-1) //channel cell below
			for(j=0;j<n_ch_per_cell;j++)
				Ap[n_circuit+nwall+i*n_ch_per_cell+j+1]++;

		if(layer!=0) //channel cell above
			for(j=0;j<n_ch_per_cell;j++)
				Ap[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j+1]++;

	}

	/*count for channel wall cells*/
	for(i=n_circuit;i<n_circuit+nwall;i++)
	{
		layer=double2int((i-n_circuit)/((n_ch_per_cell+1)*nr*nc));
		row=double2int((i-n_circuit)/((n_ch_per_cell+1)*nc))%nr;
		col=(i-n_circuit)%((n_ch_per_cell+1)*nc);

		/*regarding the silicon circuit part, the matrix has non-zero elements for silicon cells above and below the channel wall cell.*/
		col1=double2int(col/(n_ch_per_cell+1));
		Ap[layer*nr*nc+row*nc+col1+1]++; //silicon cell above
		Ap[(layer+1)*nr*nc+row*nc+col1+1]++; //silicon cell below

		/*regarding the channel wall part, the matrix has non-zero elements for itself and channel wall cells adjacent to it.*/
		Ap[i+1]++; //diagonal

		if(row!=0) //north wall cell
			Ap[i-nc*(n_ch_per_cell+1)+1]++;
		if(row!=nr-1) //south wall cell
			Ap[i+nc*(n_ch_per_cell+1)+1]++;
		if(col%(n_ch_per_cell+1)==0&&col!=0) //west wall cell
			Ap[i]++;
		if(col%(n_ch_per_cell+1)==n_ch_per_cell&&col!=nc*(n_ch_per_cell+1)-1) //east wall cell
			Ap[i+2]++;

		/*regarding the channel part, the matrix has non-zero elements for channel cells adjacent to it.*/
		col1=double2int(col/(n_ch_per_cell+1));
		index=n_circuit+nwall+layer*n_ch_per_layer*nr+row*n_ch_per_layer+col1*n_ch_per_cell+col%(n_ch_per_cell+1)-1;
		if(col%(n_ch_per_cell+1)!=0) //west channel cell
			Ap[index+1]++;

		if(col%(n_ch_per_cell+1)!=n_ch_per_cell) //east channel cell
			Ap[index+2]++;

	}

	/*count for channel part*/
	for(i=n_circuit+nwall;i<size;i++)
	{
		layer=double2int((i-n_circuit-nwall)/(nr*n_ch_per_layer));
		row=double2int((i-n_circuit-nwall)/n_ch_per_layer)%nr;
		col=(i-n_circuit-nwall)%n_ch_per_layer;

		/*regarding the silicon part, the matrix has non-zero elements for silicon cells above and below.*/
		col1=double2int(col/n_ch_per_cell);

		Ap[layer*nr*nc+row*nc+col1+1]++; //silicon cell above
		Ap[(layer+1)*nr*nc+row*nc+col1+1]++; //silicon cell below

		/*regarding the channel wall part, the matrix has non-zero elements for wall cells adjacent to it.*/
		index=n_circuit+layer*nr*nc*(n_ch_per_cell+1)+row*nc*(n_ch_per_cell+1)+col1*(n_ch_per_cell+1)+col%n_ch_per_cell;
		Ap[index+1]++; //west wall cell
		Ap[index+2]++; //east wall cell

		/*regarding the channel part, the matrix has non-zero elements for upstream channel cells and itself.*/
		Ap[i+1]++; //diagonal

		if(row!=0)
			Ap[i-n_ch_per_layer+1]++;

	}

	for(i=2;i<size+1;i++)
	{
		Ap[i]+=Ap[i-1];
	}



	printf("*************Ap Ready************\n");

	/*allocate space for Ai and Ax vector*/
	Ai=(int *)calloc(Ap[size], sizeof(int));
	Ax=(double *)calloc(Ap[size], sizeof(double));

	double rz_BEOL, rz_source, rz_bulk, rz_wall, rz_ch;
	double k_tsv=stk->k_tsv;

	double *tsv_circuitx=stk->tsv_circuitx;
	double *tsv_circuity=stk->tsv_circuity;
	double *tsv_circuitz=stk->tsv_circuitz;
	double *tsv_bondx=stk->tsv_bondx;
	double *tsv_bondy=stk->tsv_bondy;
	double *tsv_bondz=stk->tsv_bondz;

	double cw=stk->width/nr;
	double wall_width=stk->microchannel->wall_width;
	double ch_width=stk->microchannel->channel_width;

	inc=(int *)calloc(size, sizeof(int));

	for(i=0;i<size;i++)
		inc[i]=0;

	/*prepare Ai and Ax vector*/

	//Die part first
	for(i=0;i<n_circuit;i++)
	{
		layer=double2int(i/(nr*nc));
		row=double2int(i/nc)%nr;
		col=i%nc;
		//index for TSV array

		base=layer*NumSubLayers*nr*nc+row*nc+col;

		sum=0;
		/*regarding the silicon circuit part, the matrix has non-zero elements for itself and adjacent silicon cells, take note of boundary cells.*/
		Ai[Ap[i]+inc[i]]=i; //diagonal
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(col!=0) //west cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-((1+0.5*(tsv_circuitx[base]+tsv_circuitx[base-1]))*thermo->gx_BEOL[layer]
			                         +(1+0.5*(tsv_circuitx[base+nr*nc]+tsv_circuitx[base+nr*nc-1]))*thermo->gx_source[layer]
			                           +(1+0.5*(tsv_circuitx[base+2*nr*nc]+tsv_circuitx[base+2*nr*nc-1]))*thermo->gx_bulk[layer]);
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}

		if(col!=nc-1) //east cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-((1+0.5*(tsv_circuitx[base]+tsv_circuitx[base+1]))*thermo->gx_BEOL[layer]
			                         +(1+0.5*(tsv_circuitx[base+nr*nc]+tsv_circuitx[base+nr*nc+1]))*thermo->gx_source[layer]
			                         +(1+0.5*(tsv_circuitx[base+2*nr*nc]+tsv_circuitx[base+2*nr*nc+1]))*thermo->gx_bulk[layer]);
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}

		if(row!=0) //north cell
		{
			Ai[Ap[i-nc]+inc[i-nc]]=i;
			Ax[Ap[i-nc]+(inc[i-nc]++)]=-((1+0.5*(tsv_circuity[base]+tsv_circuity[base-nc]))*thermo->gy_BEOL[layer]
			                           +(1+0.5*(tsv_circuity[base+nr*nc]+tsv_circuity[base+nr*nc-nc]))*thermo->gy_source[layer]
			                           +(1+0.5*(tsv_circuity[base+2*nr*nc]+tsv_circuity[base+2*nr*nc-nc]))*thermo->gy_bulk[layer]);
			sum-=Ax[Ap[i-nc]+inc[i-nc]-1];
		}

		if(row!=nr-1) //south cell
		{
			Ai[Ap[i+nc]+inc[i+nc]]=i;
			Ax[Ap[i+nc]+(inc[i+nc]++)]=-((1+0.5*(tsv_circuity[base]+tsv_circuity[base+nc]))*thermo->gy_BEOL[layer]
			                           +(1+0.5*(tsv_circuity[base+nr*nc]+tsv_circuity[base+nr*nc+nc]))*thermo->gy_source[layer]
			                           +(1+0.5*(tsv_circuity[base+2*nr*nc]+tsv_circuity[base+2*nr*nc+nc]))*thermo->gy_bulk[layer]);
			sum-=Ax[Ap[i+nc]+inc[i+nc]-1];
		}

		/*regarding the channel wall part, the matrix has non-zero elements for wall cells above and below.*/
		if(layer!=nl-1) //wall cell below
			for(j=0;j<n_ch_per_cell+1;j++)
			{
				Ai[Ap[n_circuit+i*(n_ch_per_cell+1)+j]+inc[n_circuit+i*(n_ch_per_cell+1)+j]]=i;
				rz_source=0.5*cw/(wall_width*(1+tsv_circuitz[base+nr*nc])*thermo->gz_source[layer]);
				rz_bulk=cw/(wall_width*(1+tsv_circuitz[base+2*nr*nc])*thermo->gz_bulk[layer]);
				rz_wall=1.0/((1+tsv_bondz[i*(n_ch_per_cell+1)+j])*thermo->ch_thermal->wall_gz);

				if(j==0||j==n_ch_per_cell)
					Ax[Ap[n_circuit+i*(n_ch_per_cell+1)+j]+(inc[n_circuit+i*(n_ch_per_cell+1)+j]++)]=-0.5/(rz_source+rz_bulk+rz_wall);
				else
					Ax[Ap[n_circuit+i*(n_ch_per_cell+1)+j]+(inc[n_circuit+i*(n_ch_per_cell+1)+j]++)]=-1.0/(rz_source+rz_bulk+rz_wall);
				sum-=Ax[Ap[n_circuit+i*(n_ch_per_cell+1)+j]+inc[n_circuit+i*(n_ch_per_cell+1)+j]-1];
			}

		if(layer!=0) //wall cell above
			for(j=0;j<n_ch_per_cell+1;j++)
			{
				Ai[Ap[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]+inc[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]]=i;
				rz_BEOL=cw/(wall_width*(1+tsv_circuitz[base])*thermo->gz_BEOL[layer]);
				rz_source=0.5*cw/(wall_width*(1+tsv_circuitz[base+nr*nc])*thermo->gz_source[layer]);
				rz_wall=1.0/((1+tsv_bondz[(i-nr*nc)*(n_ch_per_cell+1)+j])*thermo->ch_thermal->wall_gz);

				if(j==0||j==n_ch_per_cell)
					Ax[Ap[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]+(inc[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]++)]=-0.5/(rz_source+rz_BEOL+rz_wall);
				else
					Ax[Ap[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]+(inc[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]++)]=-1.0/(rz_source+rz_BEOL+rz_wall);
				sum-=Ax[Ap[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]+inc[n_circuit+(i-nr*nc)*(n_ch_per_cell+1)+j]-1];
			}

		/*regarding the channel part, the matrix has non-zero elements for channel cells above and below.*/
		if(layer!=nl-1) //channel cell below
			for(j=0;j<n_ch_per_cell;j++)
			{
				Ai[Ap[n_circuit+nwall+i*n_ch_per_cell+j]+inc[n_circuit+nwall+i*n_ch_per_cell+j]]=i;
				rz_source=0.5*cw/(ch_width*(1+tsv_circuitz[base+nr*nc])*thermo->gz_source[layer]);
				rz_bulk=cw/(ch_width*(1+tsv_circuitz[base+2*nr*nc])*thermo->gz_bulk[layer]);
				rz_ch=1.0/thermo->ch_thermal->ch_gz[i*n_ch_per_cell+j];

				Ax[Ap[n_circuit+nwall+i*n_ch_per_cell+j]+(inc[n_circuit+nwall+i*n_ch_per_cell+j]++)]=-1.0/(rz_source+rz_bulk+rz_ch);
				sum-=Ax[Ap[n_circuit+nwall+i*n_ch_per_cell+j]+inc[n_circuit+nwall+i*n_ch_per_cell+j]-1];
			}

		if(layer!=0) //channel cell above
			for(j=0;j<n_ch_per_cell;j++)
			{
				Ai[Ap[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j]+inc[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j]]=i;
				rz_BEOL=cw/(ch_width*(1+tsv_circuitz[base])*thermo->gz_BEOL[layer]);
				rz_source=0.5*cw/(ch_width*(1+tsv_circuitz[base+nr*nc])*thermo->gz_source[layer]);
				rz_ch=1.0/thermo->ch_thermal->ch_gz[(i-nr*nc)*n_ch_per_cell+j];

				Ax[Ap[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j]+(inc[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j]++)]=-1.0/(rz_source+rz_BEOL+rz_ch);
				sum-=Ax[Ap[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j]+inc[n_circuit+nwall+(i-nr*nc)*n_ch_per_cell+j]-1];
			}

		Ax[diag_idx]=sum;

	}


	//Channel wall part
	for(i=n_circuit;i<n_circuit+nwall;i++)
	{
		layer=double2int((i-n_circuit)/((n_ch_per_cell+1)*nr*nc));
		row=double2int((i-n_circuit)/((n_ch_per_cell+1)*nc))%nr;
		col=(i-n_circuit)%((n_ch_per_cell+1)*nc);

		sum=0;

		/*regarding the silicon circuit part, the matrix has non-zero elements for silicon cells above and below the channel wall cell.*/
		col1=double2int(col/(n_ch_per_cell+1));
		base=layer*NumSubLayers*nr*nc+row*nc+col1;
		//silicon cell above and below
		Ai[Ap[layer*nr*nc+row*nc+col1]+inc[layer*nr*nc+row*nc+col1]]=i;
		Ai[Ap[(layer+1)*nr*nc+row*nc+col1]+inc[(layer+1)*nr*nc+row*nc+col1]]=i;

		/*circuit cell above*/
		rz_source=0.5*cw/(wall_width*(1+tsv_circuitz[base+nr*nc])*thermo->gz_source[layer]);
		rz_bulk=cw/(wall_width*(1+tsv_circuitz[base+2*nr*nc])*thermo->gz_bulk[layer]);
		rz_wall=1.0/((1+tsv_bondz[i-n_circuit])*thermo->ch_thermal->wall_gz);

		if(col%(n_ch_per_cell+1)==0||col%(n_ch_per_cell+1)==n_ch_per_cell)
		{
			Ax[Ap[layer*nr*nc+row*nc+col1]+(inc[layer*nr*nc+row*nc+col1]++)]=-0.5/(rz_source+rz_bulk+rz_wall);
			sum-=Ax[Ap[layer*nr*nc+row*nc+col1]+inc[layer*nr*nc+row*nc+col1]-1];
		}
		else
		{
			Ax[Ap[layer*nr*nc+row*nc+col1]+(inc[layer*nr*nc+row*nc+col1]++)]=-1.0/(rz_source+rz_bulk+rz_wall);
			sum-=Ax[Ap[layer*nr*nc+row*nc+col1]+inc[layer*nr*nc+row*nc+col1]-1];
		}
		/*circuit cell below*/
		rz_BEOL=cw/(wall_width*(1+tsv_circuitz[base+3*nr*nc])*thermo->gz_BEOL[layer+1]);
		rz_source=0.5*cw/(wall_width*(1+tsv_circuitz[base+4*nr*nc])*thermo->gz_source[layer+1]);
		if(col%(n_ch_per_cell+1)==0||col%(n_ch_per_cell+1)==n_ch_per_cell)
		{
			Ax[Ap[(layer+1)*nr*nc+row*nc+col1]+(inc[(layer+1)*nr*nc+row*nc+col1]++)]=-0.5/(rz_BEOL+rz_source+rz_wall);
			sum-=Ax[Ap[(layer+1)*nr*nc+row*nc+col1]+inc[(layer+1)*nr*nc+row*nc+col1]-1];
		}
		else
		{
			Ax[Ap[(layer+1)*nr*nc+row*nc+col1]+(inc[(layer+1)*nr*nc+row*nc+col1]++)]=-1.0/(rz_BEOL+rz_source+rz_wall);
			sum-=Ax[Ap[(layer+1)*nr*nc+row*nc+col1]+inc[(layer+1)*nr*nc+row*nc+col1]-1];
		}
		/*regarding the channel wall part, the matrix has non-zero elements for itself and channel wall cells adjacent to it.*/
		//diagonal
		Ai[Ap[i]+inc[i]]=i;
		diag_idx=Ap[i]+inc[i];
		Ax[Ap[i]+(inc[i]++)]=0;

		if(row!=0) //north wall cell
		{
			Ai[Ap[i-nc*(n_ch_per_cell+1)]+inc[i-nc*(n_ch_per_cell+1)]]=i;

			if(col%(n_ch_per_cell+1)==0||col%(n_ch_per_cell+1)==n_ch_per_cell)
			{
				Ax[Ap[i-nc*(n_ch_per_cell+1)]+(inc[i-nc*(n_ch_per_cell+1)]++)]=-(1+0.5*(tsv_bondy[i-n_circuit]+tsv_bondy[i-n_circuit-nc*(n_ch_per_cell+1)]))*thermo->ch_thermal->wall_gy/2;
				sum-=Ax[Ap[i-nc*(n_ch_per_cell+1)]+inc[i-nc*(n_ch_per_cell+1)]-1];
			}
			else
			{
				Ax[Ap[i-nc*(n_ch_per_cell+1)]+(inc[i-nc*(n_ch_per_cell+1)]++)]=-(1+0.5*(tsv_bondy[i-n_circuit]+tsv_bondy[i-n_circuit-nc*(n_ch_per_cell+1)]))*thermo->ch_thermal->wall_gy;
				sum-=Ax[Ap[i-nc*(n_ch_per_cell+1)]+inc[i-nc*(n_ch_per_cell+1)]-1];
			}
		}
		if(row!=nr-1) //south wall cell
		{
			Ai[Ap[i+nc*(n_ch_per_cell+1)]+inc[i+nc*(n_ch_per_cell+1)]]=i;

			if(col%(n_ch_per_cell+1)==0||col%(n_ch_per_cell+1)==n_ch_per_cell)
			{
				Ax[Ap[i+nc*(n_ch_per_cell+1)]+(inc[i+nc*(n_ch_per_cell+1)]++)]=-(1+0.5*(tsv_bondy[i-n_circuit]+tsv_bondy[i-n_circuit+nc*(n_ch_per_cell+1)]))*thermo->ch_thermal->wall_gy/2;
				sum-=Ax[Ap[i+nc*(n_ch_per_cell+1)]+inc[i+nc*(n_ch_per_cell+1)]-1];
			}
			else
			{
				Ax[Ap[i+nc*(n_ch_per_cell+1)]+(inc[i+nc*(n_ch_per_cell+1)]++)]=-(1+0.5*(tsv_bondy[i-n_circuit]+tsv_bondy[i-n_circuit+nc*(n_ch_per_cell+1)]))*thermo->ch_thermal->wall_gy;
				sum-=Ax[Ap[i+nc*(n_ch_per_cell+1)]+inc[i+nc*(n_ch_per_cell+1)]-1];
			}
		}

		if(col%(n_ch_per_cell+1)==0&&col!=0) //west wall cell
		{
			Ai[Ap[i-1]+inc[i-1]]=i;
			Ax[Ap[i-1]+(inc[i-1]++)]=-(1+0.5*(tsv_bondx[i-n_circuit]+tsv_bondx[i-n_circuit-1]))*thermo->ch_thermal->wall_gx;
			sum-=Ax[Ap[i-1]+inc[i-1]-1];
		}

		if(col%(n_ch_per_cell+1)==n_ch_per_cell&&col!=nc*(n_ch_per_cell+1)-1) //east wall cell
		{
			Ai[Ap[i+1]+inc[i+1]]=i;
			Ax[Ap[i+1]+(inc[i+1]++)]=-(1+0.5*(tsv_bondx[i-n_circuit]+tsv_bondx[i-n_circuit+1]))*thermo->ch_thermal->wall_gx;
			sum-=Ax[Ap[i+1]+inc[i+1]-1];
		}

		/*regarding the channel part, the matrix has non-zero elements for channel cells adjacent to it.*/
		col1=double2int(col/(n_ch_per_cell+1));
		index=n_circuit+nwall+layer*n_ch_per_layer*nr+row*n_ch_per_layer+col1*n_ch_per_cell+col%(n_ch_per_cell+1)-1;

		if(col%(n_ch_per_cell+1)!=0) //west channel cell
		{
			Ai[Ap[index]+inc[index]]=i;
			if(col%(n_ch_per_cell+1)==n_ch_per_cell)
				Ax[Ap[index]+(inc[index]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[index-n_circuit-nwall]+0.5/((1+tsv_bondx[i-n_circuit])*thermo->ch_thermal->wall_gx));
			else
				Ax[Ap[index]+(inc[index]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[index-n_circuit-nwall]+1.0/((1+tsv_bondx[i-n_circuit])*thermo->ch_thermal->wall_gx));
			sum-=Ax[Ap[index]+inc[index]-1];
		}

		if(col%(n_ch_per_cell+1)!=n_ch_per_cell) //east channel cell
		{
			Ai[Ap[index+1]+inc[index+1]]=i;
			if(col%(n_ch_per_cell+1)==0)
				Ax[Ap[index+1]+(inc[index+1]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[index+1-n_circuit-nwall]+0.5/((1+tsv_bondx[i-n_circuit])*thermo->ch_thermal->wall_gx));
			else
				Ax[Ap[index+1]+(inc[index+1]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[index+1-n_circuit-nwall]+1.0/((1+tsv_bondx[i-n_circuit])*thermo->ch_thermal->wall_gx));
			sum-=Ax[Ap[index+1]+inc[index+1]-1];
		}

		Ax[diag_idx]=sum;

	}


	//channel part
	for(i=n_circuit+nwall;i<size;i++)
	{

		layer=double2int((i-n_circuit-nwall)/(nr*n_ch_per_layer));
		row=double2int((i-n_circuit-nwall)/n_ch_per_layer)%nr;
		col=(i-n_circuit-nwall)%n_ch_per_layer;

		sum=0;
		ch_num=layer*n_ch_per_layer+col;
		chc_index=layer*n_ch_per_layer*nr+row*n_ch_per_layer+col;

		/*regarding the silicon part, the matrix has non-zero elements for silicon cells above and below.*/
		col1=double2int(col/n_ch_per_cell);
		base=layer*NumSubLayers*nr*nc+row*nc+col1;
		//silicon cell above
		Ai[Ap[layer*nr*nc+row*nc+col1]+inc[layer*nr*nc+row*nc+col1]]=i;
		rz_source=0.5*cw/(ch_width*(1+tsv_circuitz[base+nr*nc])*thermo->gz_source[layer]);
		rz_bulk=cw/(ch_width*(1+tsv_circuitz[base+2*nr*nc])*thermo->gz_bulk[layer]);
		rz_ch=1.0/thermo->ch_thermal->ch_gz[chc_index];

		Ax[Ap[layer*nr*nc+row*nc+col1]+(inc[layer*nr*nc+row*nc+col1]++)]=-1.0/(rz_source+rz_bulk+rz_ch);
		sum-=Ax[Ap[layer*nr*nc+row*nc+col1]+inc[layer*nr*nc+row*nc+col1]-1];

		//silicon cell below
		Ai[Ap[(layer+1)*nr*nc+row*nc+col1]+inc[(layer+1)*nr*nc+row*nc+col1]]=i;
		rz_BEOL=cw/(ch_width*(1+tsv_circuitz[base+3*nr*nc])*thermo->gz_BEOL[layer+1]);
		rz_source=0.5*cw/(ch_width*(1+tsv_circuitz[base+4*nr*nc])*thermo->gz_source[layer+1]);

		Ax[Ap[(layer+1)*nr*nc+row*nc+col1]+(inc[(layer+1)*nr*nc+row*nc+col1]++)]=-1.0/(rz_BEOL+rz_source+rz_ch);
		sum-=Ax[Ap[(layer+1)*nr*nc+row*nc+col1]+inc[(layer+1)*nr*nc+row*nc+col1]-1];

		/*regarding the channel wall part, the matrix has non-zero elements for wall cells adjacent to it.*/
		index=n_circuit+layer*nr*nc*(n_ch_per_cell+1)+row*nc*(n_ch_per_cell+1)+col1*(n_ch_per_cell+1)+col%n_ch_per_cell;

		//west wall cell
		Ai[Ap[index]+inc[index]]=i;
		if(col%n_ch_per_cell==0)
			Ax[Ap[index]+(inc[index]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[chc_index]+0.5/((1+tsv_bondx[index-n_circuit])*thermo->ch_thermal->wall_gx));
		else
			Ax[Ap[index]+(inc[index]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[chc_index]+1.0/((1+tsv_bondx[index-n_circuit])*thermo->ch_thermal->wall_gx));
		sum-=Ax[Ap[index]+inc[index]-1];

		//east wall cell
		Ai[Ap[index+1]+inc[index+1]]=i;
		if(col%n_ch_per_cell==n_ch_per_cell-1)
			Ax[Ap[index+1]+(inc[index+1]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[chc_index]+0.5/((1+tsv_bondx[index+1-n_circuit])*thermo->ch_thermal->wall_gx));
		else
			Ax[Ap[index+1]+(inc[index+1]++)]=-1.0/(1.0/thermo->ch_thermal->ch_gx[chc_index]+1.0/((1+tsv_bondx[index+1-n_circuit])*thermo->ch_thermal->wall_gx));
		sum-=Ax[Ap[index+1]+inc[index+1]-1];

		/*regarding the channel part, the matrix has non-zero elements for upstream channel cells and itself.*/
		//diagonal
		Ai[Ap[i]+inc[i]]=i;
		Ax[Ap[i]+(inc[i]++)]=sum+thermo->ch_thermal->flowrate[ch_num]*SP_FLUID;

		if(row!=0)
		{
			Ai[Ap[i-n_ch_per_layer]+inc[i-n_ch_per_layer]]=i;
			Ax[Ap[i-n_ch_per_layer]+(inc[i-n_ch_per_layer]++)]=-thermo->ch_thermal->flowrate[ch_num]*SP_FLUID;
		}

	}

	printf("*************Ai and Ax Ready************\n");



	/*initialize power matrix*/

	p=(double *)calloc(size, sizeof(double));

	for(i=0;i<n_circuit;i++)
	{
		p[i]=thermo->power_grid[i];
	}

	for(i=n_circuit;i<n_circuit+nwall;i++)
	{
		p[i]=0;
	}

	for(i=n_circuit+nwall;i<size;i++)
	{
		layer=double2int((i-n_circuit-nwall)/(nr*n_ch_per_layer));
		row=double2int((i-n_circuit-nwall)/n_ch_per_layer)%nr;
		col=(i-n_circuit-nwall)%n_ch_per_layer;

		ch_num=layer*n_ch_per_layer+col;

		p[i]=0;
		if(row==0)
			p[i]+=AMBIENT*thermo->ch_thermal->flowrate[ch_num]*SP_FLUID;


	}

	printf("*************p Ready************\n");

//	print_klu_matrix(Ap, Ai, Ax, p, size);

	klu_symbolic *Symbolic ;
    klu_numeric *Numeric ;
    klu_common Common ;

   // printf("*************KLU variable declared************\n");
    klu_defaults (&Common) ;
	Symbolic = klu_analyze (size, Ap, Ai, &Common) ;
	//printf("*************Symbolic Ready************\n");
    Numeric = klu_factor (Ap, Ai, Ax, Symbolic, &Common) ;
	//printf("*************Numeric Ready************\n");
   	klu_solve (Symbolic, Numeric, size, 1, p, &Common) ;
	printf("*************KLU solved************\n");
    klu_free_symbolic (&Symbolic, &Common) ;
    klu_free_numeric (&Numeric, &Common) ;

    /*output matrix elements to a file, for debug purpose*/
//  dump_sparse_matrix(Ap, Ai, Ax, size, "matrixdetail");


	free(Ap);
	free(Ai);
	free(Ax);

	for(i=0;i<n_circuit;i++)
	{
		thermo->steady_temp_grid[i]=p[i];
	}


	for(i=n_circuit+nwall;i<size;i++)
	{
		thermo->ch_thermal->Tch[i-n_circuit-nwall]=p[i];
	}


	/*map back to n units*/
	xlate_temp_g2b(stk,thermo);


	free(p);

}



/*map block power/temperature to grid*/
void xlate_vector_b2g(stk_descr_t *stk, thermal_data_t *thermo, double *b, int type)
{
	int i, j, n, base = 0;
	double area;
	int rows=stk->rows;
	int cols=stk->cols;
	/* area of a single grid cell	*/
	area = (stk->width * stk->length) / (stk->cols * stk->rows);

	for(n=0; n < stk->n_layers; n++)
	{
		for(i=0; i < stk->rows; i++)
			for(j=0; j < stk->cols; j++)
			{
				/* for each grid cell, the power density / temperature are
				 * the average of the power densities / temperatures of the
				 * blocks in it weighted by their occupancies
				 */
				/* convert power density to power	*/
				if (type == V_POWER)
					thermo->power_grid[n*rows*cols+i*cols+j] = blist_avg(stk->layers[n].b2gmap[i][j],
															 stk->layers[n].flp, &b[base], type) * area;
				/* no conversion necessary for temperature	*/
				else if (type == V_TEMP)
					thermo->steady_temp_grid[n*rows*cols+i*cols+j] = blist_avg(stk->layers[n].b2gmap[i][j],
										 stk->layers[n].flp, &b[base], type);
				else
					fatal("unknown vector type\n");


			}
		/* keep track of the beginning address of this layer in the
		 * block power vector
		 */
		base += stk->layers[n].flp->n_units;
	}

}

/* translate temperature between grid and block vectors	*/
void xlate_temp_g2b(stk_descr_t *stk, thermal_data_t *thermo)
{
	int i, j, n, u, base = 0, count;
	int i1, j1, i2, j2, ci1, cj1, ci2, cj2;
	double min, max, avg;
	double *b=thermo->steady_temp_block;
	double *g=thermo->steady_temp_grid;
	int rows=stk->rows;
	int cols=stk->cols;

	for(n=0; n < stk->n_layers; n++) {
		for(u=0; u < stk->layers[n].flp->n_units; u++) {
			/* extent of this unit in grid cell units	*/
			i1 = stk->layers[n].g2bmap[u].i1;
			j1 = stk->layers[n].g2bmap[u].j1;
			i2 = stk->layers[n].g2bmap[u].i2;
			j2 = stk->layers[n].g2bmap[u].j2;

			/* map the center grid cell's temperature to the block	*/
			if (thermo->map_mode == GRID_CENTER) {
				/* center co-ordinates	*/
				ci1 = (i1 + i2) / 2;
				cj1 = (j1 + j2) / 2;
				/* in case of even no. of cells, center
				 * is the average of two central cells
				 */
				/* ci2 = ci1-1 when even, ci1 otherwise	*/
				ci2 = ci1 - !((i2-i1) % 2);
				/* cj2 = cj1-1 when even, cj1 otherwise	*/
				cj2 = cj1 - !((j2-j1) % 2);

				b[base+u] = (g[n*rows*cols+ci1*cols+cj1] + g[n*rows*cols+ci2*cols+cj1] +
						g[n*rows*cols+ci1*cols+cj2] + g[n*rows*cols+ci2*cols+cj2]) / 4;
				continue;
			}

			/* find the min/max/avg temperatures of the
			 * grid cells in this block
			 */
			avg = 0.0;
			count = 0;
			min = max = g[n*rows*cols+i1*cols+j1];
			for(i=i1; i < i2; i++)
				for(j=j1; j < j2; j++) {
					avg += g[n*rows*cols+i*cols+j];
					if (g[n*rows*cols+i*cols+j] < min)
						min = g[n*rows*cols+i*cols+j];
					if (g[n*rows*cols+i*cols+j] > max)
						max = g[n*rows*cols+i*cols+j];
					count++;
				}

			/* map to output accordingly	*/
			switch (thermo->map_mode)
			{
				case GRID_AVG:
					b[base+u] = avg / count;
					break;
				case GRID_MIN:
					b[base+u] = min;
					break;
				case GRID_MAX:
					b[base+u] = max;
					break;
				/* taken care of already	*/
				case GRID_CENTER:
					break;
				default:
					fatal("unknown mapping mode\n");
					break;
			}
		}
		/* keep track of the beginning address of this layer in the
		 * block power vector
		 */
		base += stk->layers[n].flp->n_units;
	}

}

/*pad in the default flowrate*/
void fill_default_flowrate(stk_descr_t *stk, thermal_data_t *thermo)
{
	int i, total_ch_num;
	double flow_per_ch;

	total_ch_num=stk->microchannel->n_channels_per_layer*(stk->n_layers-1);
	printf("flowrate: %lf\n", stk->microchannel->default_flowrate);
	flow_per_ch=stk->microchannel->default_flowrate/total_ch_num*1e-6/60;

	printf("flow per ch is %e, default is %lf\n", flow_per_ch,stk->microchannel->default_flowrate);
	for(i=0;i<total_ch_num;i++)
	{
		thermo->ch_thermal->flowrate[i]=flow_per_ch;
	}
}


void dump_tmap(IO_t *IO, stk_descr_t *stk, thermal_data_t *thermo, char *tmap)
{
	int size=stk->rows*stk->cols*stk->n_layers;
	int layer_size=stk->rows*stk->cols;
	int i, max_index, min_index;
	double max=thermo->steady_temp_grid[0]-273.15, min=thermo->steady_temp_grid[0]-273.15;
	int n_ch_per_cell=stk->microchannel->n_channels_per_layer/stk->cols;

	for(i=0;i<stk->n_layers*layer_size;i++)
	{
		if(max<thermo->steady_temp_grid[i]-273.15)
		{
			max=thermo->steady_temp_grid[i]-273.15;
			max_index=i;
		}

		if(min>thermo->steady_temp_grid[i]-273.15)
		{
			min=thermo->steady_temp_grid[i]-273.15;
			min_index=i;
		}
	}
	printf("The highest unit temperature is %f Celcius, %f K\n", max, max+273.15);
	printf("The lowest unit temperature is %f Celcius, %f K\n", min, min+273.15);

	max_index=max_cell_temp(stk, thermo);
	printf("Overall highest temp is %f K\n", thermo->steady_temp_grid[max_index]);

	FILE *grid_temp=fopen(tmap, "wb");
	for(i=0;i<size;i++)
		fprintf(grid_temp, "%d\t%f\n", i%layer_size, thermo->steady_temp_grid[i]);
	fclose(grid_temp);

}

/*output simulation results*/
void dump_outputs(IO_t *IO, stk_descr_t *stk, thermal_data_t *thermo)
{
	int size=stk->rows*stk->cols*stk->n_layers;
	int layer_size=stk->rows*stk->cols;
	int i, max_index, min_index;
	double max=thermo->steady_temp_grid[0]-273.15, min=thermo->steady_temp_grid[0]-273.15;
	int n_ch_per_cell=stk->microchannel->n_channels_per_layer/stk->cols;

	for(i=0;i<stk->rows*stk->cols;i++)
	{
		if(max<thermo->steady_temp_grid[i]-273.15)
		{
			max=thermo->steady_temp_grid[i]-273.15;
			max_index=i;
		}

		if(min>thermo->steady_temp_grid[i]-273.15)
		{
			min=thermo->steady_temp_grid[i]-273.15;
			min_index=i;
		}

	}

	printf("The highest unit temperature is %f Celcius, %f K\n", max, max+273.15);
	printf("The lowest unit temperature is %f Celcius, %f K\n", min, min+273.15);

	if(IO->thermal_map!=NULL)
	{
		FILE *grid_temp=fopen(IO->thermal_map, "wb");
		for(i=0;i<size;i++)
			fprintf(grid_temp, "%d\t%f\n", i%layer_size, thermo->steady_temp_grid[i]);
		fclose(grid_temp);
	}

	FILE *Ch_Temp=fopen("Ch_Temp", "wb");
	for(i=0;i<stk->microchannel->n_channels_per_layer*(stk->n_layers-1)*stk->rows;i++)
		fprintf(Ch_Temp, "%d\t%f\n", i%(stk->microchannel->n_channels_per_layer*stk->rows), thermo->ch_thermal->Tch[i]);
	fclose(Ch_Temp);
}

/*clean data structures after simulation*/

void clean_thermal_data(stk_descr_t *stk, thermal_data_t *thermo)
{
	free(thermo->power_grid);

	free(thermo->gx_BEOL);
	free(thermo->gy_BEOL);
	free(thermo->gz_BEOL);
	free(thermo->gx_source);
	free(thermo->gy_source);
	free(thermo->gz_source);
	free(thermo->gx_bulk);
	free(thermo->gy_bulk);
	free(thermo->gz_bulk);

	if(stk->cooling_method==1)
		clean_pk_thermal(thermo->pk_thermal);
	else if(stk->cooling_method==2)
		clean_ch_thermal(thermo->ch_thermal);

	free(thermo->steady_temp_grid);
	free(thermo->steady_temp_block);
	free(thermo);
}

void clean_pk_thermal(package_thermal_t *pk_thermal)
{
	free(pk_thermal);
}

void clean_ch_thermal(channel_thermal_t *ch_thermal)
{
	free(ch_thermal->h_ch);
	free(ch_thermal->flowrate);
	free(ch_thermal->Tch);
	free(ch_thermal->ch_gx);
	free(ch_thermal->ch_gz);
	free(ch_thermal);
}

void clean_stk_description(stk_descr_t *stk)
{
	if(stk->cooling_method==1)
		free(stk->heatsink);
	else if(stk->cooling_method==2)
		free(stk->microchannel);

	clean_layers(stk);

	free(stk->tsv_physi);
	free(stk->tsv_circuitx);
	free(stk->tsv_circuity);
	free(stk->tsv_circuitz);
	free(stk->tsv_bondx);
	free(stk->tsv_bondy);
	free(stk->tsv_bondz);
}

void clean_layers(stk_descr_t *stk)
{
	int i,j,k;
	blist_t *ptr, *temp;
	for(i=0;i<stk->n_layers;i++)
	{
		free(stk->layers[i].flp->units);
		free(stk->layers[i].flp);
		for(j=0;j<stk->rows;j++)
			for(k=0;k<stk->cols;k++)
			{
				ptr = stk->layers[i].b2gmap[j][k];
				while(ptr) {
					temp = ptr->next;
					free(ptr);
					ptr = temp;
				}
			}
		free(stk->layers[i].b2gmap[0]);
		free(stk->layers[i].b2gmap);
		free(stk->layers[i].g2bmap);
	}

	free(stk->layers);
}

int max_cell_temp(stk_descr_t *stk, thermal_data_t *thermo)
{
	double max=0;
	int i, max_index;
	int nl = stk->n_layers;
	int nr = stk->rows;
	int nc = stk->cols;

	for(i=0;i<nl*nr*nc;i++)
	{
		if(max<thermo->steady_temp_grid[i])
		{
			max=thermo->steady_temp_grid[i];
			max_index=i;
		}
	}

	return max_index;
}

void uniform_limit_control(stk_descr_t *stk, thermal_data_t *thermo, double temp_limit)
{
	double max_temp, scale, total_flow=0;
	int i, max_index;

	max_index=max_cell_temp(stk,thermo);
	max_temp=thermo->steady_temp_grid[max_index];

	do
	{
		scale=MIN(1-6e-4*(temp_limit-max_temp)*(temp_limit-max_temp),0.97);

		for(i=0;i<stk->microchannel->n_channels_per_layer*(stk->n_layers-1);i++)
			thermo->ch_thermal->flowrate[i]*=scale;
		update_thermal_data(stk, thermo);
		solve_steady_state(stk, thermo);
		max_index=max_cell_temp(stk,thermo);
		max_temp=thermo->steady_temp_grid[max_index];

		printf("scale is %f, max_temp is %f K\n", scale, max_temp);

	}while(max_temp<temp_limit-1);

	for(i=0;i<stk->microchannel->n_channels_per_layer*(stk->n_layers-1);i++)
		total_flow+=thermo->ch_thermal->flowrate[i];

	printf("total_flow by uniform control is %f ml/min, max_temp is %f\n", total_flow*60e6, max_temp);
}

void find_max_thermal_gradient(stk_descr_t *stk, thermal_data_t *thermo)
{
	int nl = stk->n_layers;
	int nr = stk->rows;
	int nc = stk->cols;
	int i,j, row,col;
	int size=nl*nr*nc;
	double ch=stk->length/nr;
	double max_grad=0;
	int max_index1, max_index2;

	for(i=0;i<size;i++)
	{
		row=double2int(i/nc)%nr;
		col=i%nc;


		if(col!=0) //west
			if(max_grad<thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i-1])
			{
				max_grad=thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i-1];
				max_index1=i;
				max_index2=i-1;
			}

		if(col!=nc-1) //east
			if(max_grad<thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i+1])
			{
				max_grad=thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i+1];
				max_index1=i;
				max_index2=i+1;
			}

		if(row!=0) //north
			if(max_grad<thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i-nc])
			{
				max_grad=thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i-nc];
				max_index1=i;
				max_index2=i-nc;
			}

		if(row!=nr-1) //south
			if(max_grad<thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i+nc])
			{
				max_grad=thermo->steady_temp_grid[i]-thermo->steady_temp_grid[i+nc];
				max_index1=i;
				max_index2=i+nc;
			}


	}

	max_grad/=(ch*1e3);

	printf("the max_grad is %f K/mm, index1 is %d, index2 is %d\n", max_grad, max_index1, max_index2);
}

void find_max_min_power(stk_descr_t *stk, thermal_data_t *thermo)
{
	int i;
	int size=stk->n_layers*stk->rows*stk->cols;
	double max, min;
	max=thermo->power_grid[0];
	min=thermo->power_grid[0];

	for(i=1;i<size;i++)
	{
		if(max<thermo->power_grid[i])
			max=thermo->power_grid[i];
		if(min>thermo->power_grid[i])
			min=thermo->power_grid[i];
	}

	printf("the max power per cell is %f, min is %f\n", max, min);
}

void export_flowrate(stk_descr_t *stk, thermal_data_t *thermo, char *flowfile)
{
	int size=stk->microchannel->n_channels_per_layer*(stk->n_layers-1);
	FILE *fp=fopen(flowfile, "wb");
	int i;
	for(i=0;i<size;i++)
		fprintf(fp, "%e\n", thermo->ch_thermal->flowrate[i]);
	fclose(fp);
}

void export_powermap(stk_descr_t *stk, thermal_data_t *thermo, char *powermap)
{
	int size=stk->rows*stk->cols;
	FILE *fp=fopen(powermap, "wb");
	int i;
	for(i=0;i<size;i++)
		fprintf(fp, "%d\t %f\n", i, thermo->power_grid[i]*100);
	fclose(fp);
}

void load_flow_setting(stk_descr_t *stk, thermal_data_t *thermo, char *flowfile)
{
	FILE *fp=fopen(flowfile, "rb");
	int size=stk->microchannel->n_channels_per_layer*(stk->n_layers-1);
	char *ptr, line[LINE_SIZE];
	int count;
	double dval;

	fseek(fp, 0, SEEK_SET);
	count=0;
	while(count<size)
	{
		fgets(line, LINE_SIZE, fp);
		if (feof(fp))
		{
			fatal("Please Check flow file format.");
		}

		/* ignore comments and empty lines	*/
		ptr = strtok(line, " \r\t\n");
		if (!ptr || ptr[0] == '#')
			continue;

		if (sscanf(ptr, "%lf", &dval) != 1)
			fatal("please check flow rate\n");
		thermo->ch_thermal->flowrate[count]=dval;
		count++;
	}
	fclose(fp);
}
